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In ab initio calculations of electronic structures, total energies, and forces, it is convenient and often 
even necessary to employ a broadening of the occupation numbers. If done carefully, this improves 
the accuracy of the calculated electron densities and total energies and stabilizes the convergence of 
the iterative approach towards self-consistency. However, such a boardening may lead to an error 
in the calculation of the forces. Accurate forces are needed for an efficient geometry optimization 
of polyatomic systems and for ab initio molecular dynamics (MD) calculations. The relevance of 
this error and possible ways to correct it will be discussed in this paper. The first approach is 
computationally very simple and in fact exact for small MD time steps. This is demonstrated for 
the example of the vibration of a carbon dimer and for the relaxation of the top layer of the (111)- 
surfaces of aluminium and platinum. The second, more general, scheme employs linear-response 
theory and is applied to the calculation of the surface relaxation of Al(lll). We will show that 
the quadratic dependence of the forces on the broadening width enables an efficient extrapolation 
to the correct result. Finally the results of these correction methods will be compared to the forces 
obtained by using the smearing scheme, which has been proposed by Methfessel and Paxton. 
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In ab initio electronic structure and total energy calcula- 
tions the integrals over the Brillouin zone are ocmimonly 
replaced by the sum over a mesh of k-pointsEra. This 
approach is very efficient for insulators, but for metallic 
systems convergence with respect to the number of k- 
points becomes slow. Here the introduction of fractional 
occupation numbers is a convenient way to improve the 
k-space integration and in addition to stabilize the con- 
vergence in the iterative approach to self-consistency. In 
these broadening schemes the eigenstates are occupied 
according to aj^maoth function, e.g. a Gaussiantl or the 
Fermi functionQu'LTa at a finite temperature. 
When a broadening scheme is employed in a density func- 
tional theory calculation, the computed electron density 
of the ground state n.o(r) does not minimize the func- 
tional of the total energy E but the functional of the free 
energy A: 



A [Tci; n]^ E [Tci; n] - T^S [T^i; n] 



(1) 



where S denotes the entropy associated with the occupa- 
tion numbers of the Kohn-Sham orbitals and Td is the. 
broadening parameter. In the case of Fermi broadeningEl 
we get: 



s^-kBj2 [/' 111 + (1 - In (1 - m 



(2) 



Since the temperatures commonly used for the broaden- 
ing are much higher than the physical ones (it is con- 
vinient to use fesTd ^ 0.1 eV), neither the total energies 
nor the free energies (Eq. |l]) are directly meaningful. 
One way to obtain the ground state energy at zero tem- 
perature is based on the well known factS, that for the 



free electron gas, the quantities A and E depend quadrat- 
ically on Td. Therefore one can write: 

A{Tci)^E^°^°-Ij'T^, + 0{T^,) (3) 



^zcro 
^zcro ^ 



E (Tel) - + ^I'T^i + O (r<?i) . (4) 

As it was pointed out by GillanB, it follows from Eqs. (|^), 
(^ and (^) that the extrapolation of the total energy 
towards the Tei=0 result is: 



^EiTci)-^TciS{Tci) . 



(5) 
(6) 



Obtaining E^'"-'" using Eqs. (^ and (|^) is straightforward 
and gives very satisfactory resultsu. This is shown in Fig- 
ure H(a) for a slab consisting of four layers of Aluminium. 
Here the extrapolation (filled circles) matches perfectly 
the zero-temperature energy, even for quite large broad- 
ening temperatures. For a systems like Platinum, which 
was chosen as an example for a not free electron like 
system. Figure |l](b) shows, that this extrapolation is in- 
deed an approximation. But if the broadening parameter 
is chosen carefully (typically used broadenig parameters 
are about 0.1 eV or lower) the extrapolation still gives 
acceptable results. 

Calculating the forces, however, is more complicated. 
The forces on atoms are defined as the derivative of the 
total energy E^'''^° with respect to the atomic positions: 

^1 rpzcvo 

pzcio ^ 



dA 1 dS 



(7) 
(8) 
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FIG. 1. Total energy ii'(fcBrci) (open diamonds, 

dash-dotted line) and energy extrapolated to fceTci = eV 
^zcro ^gjjg(j circles, solid line) as a function of k^T using a 
Fermi broadening and total energy (open squares, dashed 
line) using the Methfessel-Paxon smearing of first order for a 
slab of four-layers of aluminuim (a) and platinum (b). Dots 
present computed values, lines are fits to guide the eye. 



While the quantity 



dA 



(9) 



is easily evaluated, the evaluation of dS/dR is some- 
what more elaborate. Neglecting the entropy term in 
Eq. (^) implies that the forces are not in agreement 
with the gradient of the total energy of Eq. (||). As 
a consequence, when those forces are used to relax 
the atoms towards their equilibrium positions, the ob- 
tained geometry is likely to be different from that which 
minimizes E^°'^°. Figure ^ presents results for a four- 
layers aluminum (111) slab,p5idth fully separable, norm- 
conserving pseudopotentialsOr't, .a plane-wave basis set 
(£"="* = 8 Ry) and 18 k-pointsllil to sample the surface 
Brillouin zone. An untipically high broadening of 0.5 eV 
was used to show the effect and it is clearly visible that 
the Hellmann-Feynmann- Force F^^(Tci) acting on the 
surface layer vanishes for the position minimizing the free 
energy, as it should according to Eq. (|^). But the mini- 
mum of the total energy is at a different position. 
Even at this untypically high broadening temperature, 
which was chosen to illustrate the effect, the error in the 
equilibrium position is only 10~^ nm, which is less than 
one percent of the interlayer distance. This indicates, 
that the error in the forces due to occupation number 
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FIG. 2. (a) Zero-temperature energy and (b) free en- 

ergy A of a four-layers Al (111) slab as a function of the 
height Z of the first atomic layer. To demonstrate the effect, 
a Fermi-Dirac broadening with an untypically high broaden- 
ing pararmeter fcsTei^O.S eV has been used, (c) Force acting 
on the first atomic layer (open dots are non-corrected and 
full dots are corrected forces). Dots present computed values, 
lines are fits to guide the eye. 

broadening may be neglegible in the case of relaxations 
if the temperature is chosen not unreasonable high, but 
unfortunally this has to be checked for each particular 
system. 

The situation is less clear and the needs are more de- 
manding, when it comes to molecular dynamics (MD) 
simulations. When the non-corrected forces (Eq. 0) are 
used to perform MD, the sum of potential kinetic energies 
of the atoms is not a conserved quantity. An illustration 
of this is shown in Figure ||(a) for a MD simulation of the 
vibration of a carbon dimer. The graph shows the total 
energy (open circles) and the free energy (closed circles) 
of the systems as a function of the interatomic distance. 
The trajectory was integrated over approx. two periods 
using the verlet algorithm and no thermostat had been 
employed. Using the non-corrected forces produces a mo- 
tion in which the total energy is not conserved but oscil- 
lates with the frequency of the motion. Figure |^(a) cleary 
shows, that at the turning points, where the kinetic en- 
ergy vanishes, the potential energies are different, which 
is obviously unphysical. In fact, the free energies, given 
by Eq. (^, at the turning points are equal. 
Plotting the eigenvalues as a function of the intermolec- 
ular distance (Fig. 0(b)) shows a crossing of the 2a g level 
and the two-fold degenerated l7r„ levels near 
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FIG. 3. Free Energy E''"^° (filled circles) and total energy 
E (open circles) as a function of interatomic distance during 
the vibration of a carbon dimer. The points are taken at dif- 
ferent times during the vibration. The border points at the 
outer left and right are the turning points of the motion, thus 
the range of distances covered in each figure correspond to 
the amplitude of the vibration, (a) Non-corrected forces have 
been used for the molecular dynamics. Obviously the poten- 
tial energy at the two turning points is very different, (b) 
Corrected forces have been used for the molecular dynamics. 
The potential energy at the two turning points is nearly the 
same. 

the equilibrium distance do- Since in the groundstate, 
the 2(Tg level is empty (LUMO), while the ItTu levels are 
fully occupied (HOMO), the introduction of a broaden- 
ing function causes a noticeable change in the occupation 
numbers of these orbitals (Fig. ^a)) during the vibra- 
tion. This leads to a non- negligible correction to r(rci) 
according to Eq. ||. We will now describe how this en- 
tropy contribution to the forces can be evaluated. From 
Eq. (^, one obtains the expression for TddS/dR: 



dS_ 
'9R 



an 1 - /, 



(10) 



which, in the case of a Fermi distribution for the occu- 
pation numbers, reduces to 



dS_ „ g/» 



OK 



(11) 



The ei denote the energies of the Kohn-Sham orbitals. 
Using a Fermi distribution for the fi, we find, that in the 
case of the vibrating dimer on the relevant length scales 



the occupation numbers change linearly with R (see 
Fig. 11(a)), i.e. dfi/dli is nearly a constant. Figure ^(b) 
illustrates the symmetric and total energy-conserving vi- 
brations obtained when the forces are corrected according 
to Eqs. (^) and (|ll|) using the linear treatment, which is 
well satisfied for the small time steps which are typically 
used in a MD simulation. We expect that for many ab 
initio MD calculations a good estimate of dfi /dH can be 
obtained from the past atomic geometries. This is partic- 
ularly true because due to the small time steps {At ~ ^ 
of the period of the vibration) the difference between ad- 
jacent geometries in a MD calculation is typically very 
small, i.e. only about 10"'* nm. 
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FIG. 4. (a) Occupation numbers of the highest occupied 
molecular orbital lag (HOMO, two- fold degenerate) and the 
lowest unoccupied molecular orbital (LUMO) 2-Ku as a func- 
tion of the interatomic distance during the vibration of a car- 
bon dimer, using fcBTei=0.5 eV. 

(b) Energy levels of a carbon dimer as a function of the in- 
teratomic distance, do marks the quilibrium distance. 

To complete our analysis, we also applied a more general 
method, for calculating dfi/dH. We will use again the Al 
slab as an example. At constant temperature and number 
of electrons, the occupation numbers fi depend on the 
one-electron energies and on the chemical potential fi. 
Thus, 



dfi _ dfi / dej dfi 



(12) 



in which the chemical potential is obtained from the con- 
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straint that the fi sum up to the number of electrons. 
Because a given atomic displacement will result in an in- 
crease of some eigenvalues and a decrease of the other 
ones, the derivative of the chemical potential is smaller 
than those of the eigenvalues. It is nevertheless not neg- 
ligible: in the case of the relaxation of Al (111), we find 
that the contribution of the derivative of the chemical 
potential is about 0.2 times that of the derivative of the 
eigenvalues. _ _ 

Linear-response theoryEjlij enables us to calculate the 
quantities dei/dTl, which are the expectation values of 
dH/dR for the eigenstates \4>i) 



OH 



(13) 



Other methods would be more time-consuming and thus 
inadequate for the purpose of geometry optimization or 
MD, where the forces must be calculated for many dif- 
ferent atomic configurations and sometimes for systems 
of 100 atoms and more. The dielectric matrix, for exam- 
ple, would require the calculation of many bands and the 
inversion of large matrices. 

The explicit dependence of Ti. on R is only in V^xt, but Vh 
and Vxc depend on the atomic positions through the elec- 
tron density. For that reason, an evaluation of the forces 
not using the values of the of previous atomic positions 
requires a self-consistent calculation of the derivatives 
dei/dR. Neglecting the dependence of Vh and V^c on R 
to avoid the use of an iterative computation would result 
in unacceptable errors: in the example of the Al slab, we 
encountered cases where non self-consistent forces were 
too large by an order of magnitude. 
The derivative of the electron density is given by 



dn{r) 
dR 



dR 



(14) 



where the first term accounts for the redistribution of the 
electrons among the orbitals due to the variation of the 
ei and the second term comes from the modification of 
the wavef unctions. From the current approximation of 
dei/dR, Eq. ( p^ ) is used to calculate the first term. The 
second term requires the resolution of 
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(15) 



Equation ( |l5|) is ill-conditioned, because the operator 
{Ti. — ei) of the left-hand side has in general eigenvalues 
of either sign, some of them having small absolute val- 
ues. Later in this paper (Eqs. |l8| and ^ , we are going to 
explain an iterative resolution method instrumental for 
positive-definite operators and which converges better if 
the eigenvalues of that operator are large. In order to 
make that method applicable to our problem, we sep- 
arate the Hilbert space in two subspaces: the first one 



is spanned by the computed eigenstates, and the second 
one is its complementary subspace (spanned by unoc- 
cupied states). In practical calculations, only the i,nax 
lowest energy levels, including all occupied states and 
some unoccupied states, are computed and the occupa- 
tion numbers are fractional only for the levels with ener- 
gies ztfesToi around the Fermi level. Let V be the pro- 
jector on the second subspace [i > imax) and dipi/dR 
be V Idipi/dR). Now, we can rewrite Eq. (|lj) 



as 



dn{r) 
dR 



2 i^J^di 





dH 




dR 



+2m 



^i(r) 



dR 



(16) 



Some of the differences (e^ — ej ) appearing in the denom- 
inator of the second term in Eq. (|l6|) are small but, in 
that case, the difference between the occupation numbers 
in the numerator is small as well, so the whole fraction 
has a finite value. 

The quantities dipi/dR) are the solutions of 



in - e,) 



(17) 



the 



Since Eq. (|17| ) is only for the subspace i > i„ 
operator (Ti — e^) is positive-definite. Moreover, an ap- 
proximation T) to (7i — ei) can be obtained by neglecting 
its off-diagonal matrix elements. In the plane-wave basis 
set we use in the examples, this amounts to including 
the kinetic part of the hamiltonian and the average value 
of the effective potetitial. The algorithm introduced by 
Williams and Solertj can be generalized to solve Eq. 
iteratively. An initial guess is given by 



d^i. 



(0) 



dR 



(18) 



and the sequence of approximations is obtained using 



(n+l) 



dR 



d^\ 



in) 



dR 



{e^ - H) 



d^. 



(") ' 



dR 



(19) 
(20) 



In order to make the convergence as fast as possible, 
the largest value keeping the algorithm stable is chosen 
for the constant A. The operators 'D~^ and e^^'^ are 
easily computed, because V is diagonal. The operator 



4 



is different at each iteration, because the value 
of 9n(r)/9R on which it depends is updated each time, 
in order to make the charge redistribution converge to 
self-consistence. This is analogoU|S-|to the method used 
for the total-energy minimizatior£3, where the conver- 
gence towards the self-consistent charge density and the 
diagonalization of the hamiltonian are performed simul- 
taneously. 

The method explained above was applied to the Al slab. 
The corrected force is plotted in Figure ||(c) and it in- 
deed vanishes for the geometry that minimizes the zero- 
temperature total energy. 

As expected from Eqs. (||) and the dependence of 
the non-corrected force on Tei is quadratic. Therefore it 
should be possible to extrapolate the force to the value at 
zero temperature from two points at finite temperature. 
Figure 1^ (a) shows the result for the Al(lll) slab, using 
48 k-points in the irriducible part of the Brillouin zone. 
Using the formula 



dA{Ti)/dR ~ {Ti/T2fdA{T2)/dK 



{Ti/T2f 



1 



(21) 



we obtain a value which matches quite well the nearly 
constant value, obtained for the entropy-corrected force 
and by the linear treatment. This extrapolation scheme 
even works for metals which are not free electron like, as 
it is shown in Figure ^ (b) for the case of platinum. In 
the case we used 69 k-points in the irriducible part of 
the Brillouin zone. As for a transition metal no broad- 
ening temperature above 0.1-0.2 eV should be used to 
retain the properties of the Fermi surface, in many cases 
there might be no need to correct the forces at all. But 
calculating the forces at two different temperatures and 
extrapolating to T^i = may lead to some additional 
improvement. i— ■ 

Methfessel and PaxtorJlH proposed an improved smearing 
scheme in which the orbitals are occupied according to 
a smooth approximation of the step function. The com- 
puted quantities (energy, forces,...) converge towards 
their zero-temperature values when the order of the ap- 
proximation is increased. Unfortunately higher-order ap- 
proximations become more and more wigglier and there- 
fore require larger k-point sets. In practical calculations 
therefore only the first order approximation is used. In 
Figure |l| the differences of the calculated energies for fi- 
nite temperatures to the energy at zero temperature for 
the Al-111 (a) and the Pt-111 (b) slab is plotted as a 
function of the broadening parameter. For both systems 
the deviation of the energies obtained by using the first 
order Methefessel-Paxton scheme (full circles) ist compa- 
rable to that of the extrapolated energies (open squares) 
using the finite temperature approach (open diamonds). 
Figure ^ shows that for both systems the error in the 
forces which are obtained when the Methefessel-Paxton 
scheme is used also is comparable to the error in the 
forces obtained by the quadratic extrapolation. But Fig- 
ure H also shows, that the forces obtained by using the 



MP-scheme are wigglier, especially for the Al-111 slab, 
where only 48 k-Points have been used. 
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FIG. 5. Non-corrected force acting on the first layer 
of a four layer (111) slab of aluminium (a) and plat- 
inum (b) as a function of the broadening energy fcsTei 
using a fermi-smearing (closed circles) and the Methfes- 
sel-Paxton-scheme of first order (closed squares) . Parabolic fit 
(dotted line) and forces extrapolated to T=0 (open markers) . 
Dashed line corresponds to the minimum of the parabola. 

In conclusion, we have shown how the entropy arising 
from the broadening of the occupation numbers can be 
included in the calculation of the forces. For small distor- 
tions the dependences of the occupation numbers are lin- 
ear to a good approximation. Thus, information from the 
history of the MD run can be used to determine dfi/dH. 
For the general case, e.g. when the needed information is 
not available from the history, we developed an iterative 
method based on the linear-response theory. The forces 
obtained in this way are the exact derivatives of the ex- 
trapolated zero-temperature energy. Neglecting the con- 
tribution of Vh and 14c to dei/dR, in other words stop- 
ping the calculation after the first iteration would result 
in unacceptably errors. If k-point convergence is fulfilled, 
the corrected force is nearly independent of the broaden- 
ing temperature in a wide temperature range. In that 
case, extrapolation to zero temperature from the results 
at two finite temperatures also gives good results. 
Among the three methods presented in this paper, linear- 
response theory is the most general one: it is valid for any 
dependence of the one-electron energies on the atomic 
positions. On the other hand, this method is compu- 
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tationally very expensive, which limits it's uscfulencss 
in practical calculations. The extrapolation method re- 
quires the calculation of the forces for two different tem- 
peratures, which in general does not double the number 
of iterations, since the number of iterations needed is 
much smaller starting from a selfconsistent charge den- 
sity at a different temperature. Therefore extrapolation 
from two different temperatures is more practical, espe- 
cially if the number of degrees of freedom is large. 
While the error in the forces is small compared to the 
usual accurancy in atomic relaxations, and therefore 
might be neglected if the broadening temperature is 
choosen carefully, the correction is especially important 
to stay consistent in a molecular dynamics simulation, 
where the quantity to be conserved should be the ground 
state energy of the electronic system and not the unphys- 
ical free energy of the electronic system, which is excited 
due to the broadening. 

Using the scheme proposed by Methfessel and Paxton 
in it's first order approximation leads to energies, which 
are comparable to the energies obtained by extrapolating 
the finite temperature energies to Tgi = 0. The calculated 
forces are the derivatives of these energies and need no 
further correction. 
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